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ABSTRACT 
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We  considerrthe  problem  of  computing  the  station2u:y  flows  of  a 
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viscoelastic  fluid  flowing  through  a  given  domain*  nie  proposed  numerical 
method  is  based  on  optimal  control  techniques,  which  replace  the  original 
equations  of  the  problem  by  a  minimization  problem  to  be  solved  by  a  descent 
method.  Such  techniques  are  very  powerful  and  can  handle  equations  which 
change  type,  provided  that,  as  done  here,  one  uses  2m  adequate  preconditioning 


strategy  and  that  one  computes  efficiently  the  gradient  of  the  function  to  be 
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minimized. 
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SIGNIFICANCE  AND  EXPLANATION 


Ite  consider  the  problem  of  computing  the  stationary  flows  of  a 
viscoelastic  fluid  flowing  through  a  given  domain.  The  proposed  numerical 
method  is  based  on  optimal  control  techniques#  which  replace  the  original 
equations  of  the  problem  by  a  minimization  problem  to  be  solved  by  a  descent 
method.  Such  techniques  are  very  powerful  and  can  handle  equations  which 
change  type#  provided  that#  as  done  here#  one  uses  an  adequate  preconditioning 
strategy  and  that  one  computes  efficiently  the  gradient  of  the  function  to  be 
minimized. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
svsmary  lies  with  MBC#  and  not  with  the  author  of  this  report. 


OPTIMAL  CONTROL  TECHNIQUES  FOR  COMPUTING 
STATIONARY  PLOWS  OP  VISCOELASTIC  PLUIDS 


Patrick  Lt  Tallcc* 


1.  INTRODUCTION 

The  objective  of  this  paper  is  to  present  from  a  basic  numerical 
point  of  view  a  new  class  of  methods  for  the  numerixal  calculation  of 
viscoelastic  flows.  These  methods  consist  in  : 

(i)  rewriting  the  governing  equations  as  a  least-squares  problem. 
Here,  it  is  critical  to  use  the  right  norms  and  to  introduce 
a  preconditioning  operator.  For  example,  working  with  the  quan¬ 
tities  |div  cr  -  ^  P  dx  is  completely  inadequate  and  is  wrong 
from  a  functional  analysis  point  of  view.  In  this  paper,  the 
preconditioning  strategy  will  rely  on  the  introduction  of  a  pi¬ 
vot  space; 

(ii)  solving  the  resulting  minimization  problem  by  a  Finite  Element 
Method  associated  to  a  conjugate  gradient  or  to  a  Quasi-Newton 
algorithm. 

This  paper  first  introduces  the  governing  equations  of  the  problem 
(§2)  and  rewrites  them  as  a  least-squares  (or  optimal  control)  problem 
(§3) .The  practical  calculation  of  the  cost  function  and  of  its  gradient 
is  then  discussed  (§4)  and  several  descent  methods  for  solving  the  resul¬ 
ting  minimization  problem  are  presented  (§5).  Some  details  on  the  practical 
implementation  of  these  ideas  and  numerical  results  are  finally  presented 
(§6  and  §7). 
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The  framework  of  the  paper  will  be  rather  general  and  could  be  used 
for  the  derivation  and  the  study  of  different  classes  of  numerical  methods. 
Moreover,  although  introduced  on  a  standard  steady  state  problem,  it  can 
be  easily  adapted  to  the  solution  of  true  evolution  problems,  of  pseudo 
evolution  problems  (time  marching  techniques)  or  to  an  arc  length  continua¬ 
tion  approach  (GLOWINSKI  [  1984,  p.206  ]).  In  any  case,  in  the  proposed  siethods, 
approximation  errors  can  be  made  very  small  and  change  of  type  in  the  gover¬ 
ning  equations  should  not  be  damaging  as  seen  from  the  analogy  with  transonic 
flow  computations. 


2.  EQUATIONS  OF  THE  PROBLEM 

Let  us  consider  a  viscoelastic  fluid  (say  of  upper-convected  Maxwell 
type),  flowing  viscously  through  a  given  domain  (Fig.1).  We  suppose  that 
no  slip  occurs  along  the  solid  walls  and  that  the  fluid  velocity  at  the 
entrance  and  at  the  exit  of  the  domain  is  given. 


Pr 

Figure  1  -  The  physical  problem 
(out  of  Malkus  (19841) 
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For  example,  such  situations  arise  while  studying  plane  flows  over 
slots,  such  as  those  studied  experimentally  by  Bird  and  al.  [  1982  ] .  The 
equations  governing  such  situations  are  simply  : 

( I )  EQUILIBRIUM 


-div(o)  +  p(u*7)u  “  f  in  H  , 

(2)  CONSTITUTIVE  LAW  (upper-convected  Maxwell) 


£  ■  2d  ■  Pi  • 

0jj(x,t)  -  /"  • 

-o»  X  ^ 

<  3X,.(x,r) 


Xj(x,T)  ■  position  at  r  of  the  particle  which  is  in  x  at  time 
t  and  which  is  subjected  to  the  velocity  field  u. 


(3)  KINEMATIC  RESTRICTIONS 

div  u  •  0,  u  ■  Uj  on  r  . 

Above,  u  represents  the  fluid  velocity,  o  the  Cauchy  stress  tensor,  p  a 
hydrostatic  pressure,  p  the  fluid  density,  y  the  fluid  viscosity  and  X 
the  relaxation  time.  Observe  that,  as  an  extra  boundary  condition,  the 
constitutive  law  (2)  requires  the  knowledge  of  what  happened  to  the  fluid 
before  it  enters  the  domain. 


It  has  been  observed  in  Joseph,  Renardy,  Saut  [  1984  ]  ,  that  these 
equations  change  type  when  the  viscoelastic  Mach  number  U/Zy/Xp  reaches  1 
(U  is  a  characteristic  velocity  of  the  considered  flow).  Real  characteristics 
then  appear  along  which  the  vorticity  can  be  discontinuous.  However,  most 
numerical  methods  employed  for  solving  (l)-(3)  (such  as  the  classical  fixed 
point  method  solving  iteratively  for  the  velocities  and  then  for  the  stresses) 
cannot  handle  this  change  of  type. 
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The  idea  of  this  paper  is  to  esiploy  for  Che  nuaierical  solution  of 
(l)-(3)  optimal  control  techniques  in  a  dual  space,  which  were  used  with 
success  in  transonic  flow  computations  (GLOWINSKI  [  1984  ]),  where  a  similar 
change  of  type  occurs. 


3.  LEAST-SQOARES  FORMULATIONS  OF  THE  PROBLEM. 


3.1.  A  one-dimensional  model  problem.  Let  f  :  K  K  be  differentiable 
and  consider  the  problem  of  solving  nttmerically  the  nonlinear  equation 


f(x)  -  0  . 


If  it  has  a  solution  in  K,  then  this  equation  is  equivalent  to 


Minimise  -  |f(x)l^  over  R,  (a  >  0)  , 


problem  which  can  be  numerically  solved  by  Che  gradient  algorithm 


+  *0  -  given  , 

+  for  n  >  0  until  satisfied  set 


x  -x  f(x)f'(x) 

n+1  n  a  n  n 


This  algorithm  can  be  a  very  efficient  method  for  solving  f(x)  •  0,  provided 
Chat  a  is  properly  chosen  and  Chat  f(x)f'(x)  is  easy  to  compute.  It  will  be 
Che  basis  of  the  nusierical  technique  that  we  will  use  to  solve  (l)-(3). 


3.2.  Maxwell  viscoelasticity  in  primal  variables.  Let  V  be  the  space 
for  the  unknown  velocity.  Here,  we  Cake 


{w€  (n)  ,  div  w  -  0)  , 


whose  topological  dual  we  denote  by  V*.  Let  us  also  introduce  the  auxiliary 
unknown  z  •>  Pu,  P  being  a  given  isomorphism  from  V  onto  a  Hilbert  space  H, 

H  with  scalar  product  (.,.)  and  H  identified  to  its  dual.  The  interest  of 
Che  auxiliary  unknown  is  that,  for  P  properly  chosen,  the  viscoelastic  pro¬ 
blem  may  be  better  conditioned  with  respect  to  this  new  unknown,  and 


therefore  easier  to  solve  numerically.  In  other  words,  P  is  a  precon¬ 
ditioning  operator.  For  example,  one  can  think  of  defining  P  u  as  curl 
u,  H  as  Che  image  of  V  by  P  and  (.,.)  as  the  scalar  product.  The 
vorticicy  curl  u  will  then  be  the  pristary  variable,  which  might  be  a 
very  good  choice  because  vorticicy  is  the  canonical  variable  appearing 

in  Che  analysis  of  change  of  type.  Another  choice  of  P,  for  which  again 
T 

P  P  18  equal  to  Che  Stokes  operator,  will  be  proposed  at  Che  finite  element 
level  in  §5. 

Finally,  once  V  and  P  defined,  let  us  introduce 
L  :  V  H  ,  L(w)  -  f  .jrfx  , 

T  :  V  -*  V*, 

Vw  €  V,  u  •  Uj+  u^,  20^5^  being  given 
through  Che  constitutive  relation  (2) . 

With  these  notations,  summarized  in  the  following  diagram 


equations  (1)  to  (3)  take  the  form 

(4)  T(u  -  u,)-L  «  0  in  V*  ,  u  -  u,  €V  . 


If  (4)  has  a  solution,  then  it  is  obviously  equivalent  to  the  H  least- 
squares  formulation  : 


(5) 


MINIMIZE  J(v)  -  I  (P^Cv),  P^^v))  over  H 
WHERE  ;^(v)  £  V  IS  THE  SOLUTION  OF  THE  LINEAR  PROBLEM 
(P;^(v),  P(w))-<T(p“'v),w>-L(w),  Vwev, 


Indeed,  if  (4)  has  a  solution  u  =u-u,,  and  if  we  set  v  •  z  ~  P(u  ) 
in  (5) ,  then  the  right-hand  side  of  (5)  is  equal  to  zero,  thus  the  asso¬ 
ciated  state  vector  y(z  )  is  also  equal  to  zero  and  therefore  J(z  )  is 
equal  to  zero.  Since,  by  definition,  J(.)  is  always  positive  on  H,  this 

implies  that  z  is  a  minimizer  of  J  over  H. 

~o 

Conversely,  let  ^  be  a  minimizer  of  J  over  H.  As  seen  above,  since 
(4)  has  a  solution,  J  attains  the  value  0  on  H  and  therefore  J(z  )  must  be 
equal  to  0.  By  definition  of  J,  this  can  only  happen  if  i®  equal  to  0, 

that  is  if  T(P  'z  )-L“0  in  V*,  which  means  that  P  * (z  )  is  then  a  solution 
of  (4), 

In  summary,  if  we  assume  the  existence  of  a  solution  Co  our  original 
problem  (l)-(3),  we  can  replace  these  equations  by  the  equivalent  minimiza¬ 
tion  problem  just  written  above.  This  minimization  formulation  is  the  one 
which  will  be  used  in  our  numerical  techniques.  It  reduces  our  initial  pro¬ 
blem  to  an  optimal  control  problem,  if  we  identify  v  to  a  control  variable, 
to  a  state  vector,  (5)  to  a  state  equation  and  J(*)  to  a  cost  function. 

Here,  the  state  vector  belongs  to  the  velocity  space,  which  relates 
to  some  extend  this  formulation  to  the  one  used  by  TANNER  [  1985  ]  ,  or  by  the 
other  authors  for  which  the  velocity  is  the  working  variable. 

3.3.  Arclength  continuation  in  mixed  variables.  To  illustrate  the 
many  directions  in  which  the  above  optimal  control  framework  can  be  applied, 
we  will  briefly  and  formally  discuss  another  least  squares  formulation  of 
problem  (I)-(3).  It  is  based  on  the  approach  of  CROCHET  and  al  [  1985)  or  of 
6ERIS,  AMSTRONG  and  BROWN  [  1985  1  ,  among  others,  in  which  the  original  pro¬ 
blem  is  written  in  mixed  variables  (velocity  +  added  stress)  and  is  put  in¬ 
side  an  arclength  continuation  framework. 


Vn>  0.  find  X^)  €  V  x  I  x  H  such  that 

»  .  ■  #so»#so»  *■ 

<  >•<>.*  Is-'J  «  » -  X  . 

(P  u”  -  P  u”  fP  u”  -  P  u"  )  +  (a”  -  o*'  ,0”  “  2!!  ) 

'MO  'MO  *  'MO  'MO 

+  (X*'  -  x”  )  “  As^  , 

from  which  we  derive  our  final  lease-squares  formulation  of  the  arclength 
continuation  problem 

f  V  n>0.  Minimize  the  cost  function  (residual  dual  norm) 


I  over  the  control  space  K  x  H  x  E,  the  state  vector  (residual) 
(7)  I  in  V  X  X  X  11  being  defined  through  the  state  equation 


n-1  n-1  n-2 


n-'  n-^  n-*. 


.  /%  A  V  y  w  \  y  w  \ 

U  »  As  -(X-X  )  -(z  -z  ,z  ,z  )-(a  -a  ,a  --a  )  • 

'  '  ^  ''MO  'M) 


Observe  that  (7),  as  it  must  do,  deals  with  the  right  dual  norms  of  the 
residuals,  and  that  it  is  in  fact  equivalent  to  the  standard  discretiza¬ 
tions  of  (6)  used  in  the  litterature. 

4.  CALCULATION  OF  THE  COST  FUNCTION  AND  OF  ITS  GRADIENT 


The  minimization  formulations  of  §3,  although  equivalent  to  classical 
ones,  are  very  interesting  because  they  can  be  solved  by  a  different  class  of 
stable  numerical  algorithms.  Indeed,  any  available  descent  method  (conjugate 
gradient,  Buckley-Lenir,  Quasi-Newton...)  can  be  successfully  used  for  their 
numerical  solution.  The  practical  implementation  of  such  methods  only  requi¬ 
res  the  knowledge,  at  given  controls  v,  of  the  cost  function  J(v)  and  of  its 
gradient  J'(j^. 


Now,  the  cost  function  being  a  quadratic  function  of  the  state  vector, 
the  latter  being  the  image  of  a  nonlinear  complicate  operator  acting  on  the 
control  vector,  the  computation  of  J  or  of  J'  is  not  an  easy  matter.  Never¬ 
theless,  even  in  Che  more  complex  situation  of  §3.2,  and  unlike  a  classical 
Newton  method  which  would  require  O(N^)  operations  (N»dim  V)  to  compute  this 
gradient,  the  computation  of  J'  can  be  done  in  0(N  )  operations  by  introdu¬ 
cing  an  adjoint  state  vector  which  reduces  this  computation  to  the  explicit 
integration  of  functions  locally  defined  on  the  support  of  each  trial  func¬ 
tion. 

To  see  that,  in  the  framework  of  §3.2,  let  us  introduce  the  adjoint 
state  G(v)  defined  as 


(8)  £(v)(x,t)-i/  exp(^-^)(F^(x,r))Vx“(x,T))F^(x,t)dT 


D  -  0  if  X  e  n  , 

(5)  \  3x; 

It  •  IT  • 

X^(x,i')  ■»  position  at  time  t  of  the  particle  which 
is  in  ^  at  time  t  and  which  is  subjected 
to  the  velocity  field  ju  (x)  “  -uC*)  • 

This  adjoint  state  can  be  computed  by  an  explicit  integration  along  the 
trajectories  of  u“P  (v)+u^.  Then,  we  can  prove 

THEOREM  :  The  gradient  J ' iy)  is  the  solution  of  the  linear  problem 

(10)  ( J '  (v) ,  EV)  -  ^  p {  (w.V)u  +  (u.  V)w)  .X(v)dx 

-  ^(s-£>^>-£(v)dx,vwev. 

Proof  :  By  definition  of  the  gradient,  we  have 


(J'()^,z  )  =  lim^[J(v  +  tz)  -  J(v)l  =  (P^,  P6;^) 
~  t-O  ~  ~ 


-  .  .  ’ 

.  t  \  V.  *  t  L  *  .  •  _  «  .  w 

: /:  >>•  VrVrlry..' ’.-y  v'- 

where  ^  ■  ^(v)  is  the  solution  of  the  state  equation  (5)  and  where  8^ 
is  obtained  from  ^  by  differentiation  of  (5),  that  is  by  solving 


_1  _i 

(P6;^,  PW)  =<T'(P  v).P  z,  w>,  Vw€V  . 


Substituting  this  definition  of  in  the  expression  of  the  gradient, 
we  get 


(J'(v).z)  =<T'(P-*V).P-*Z,  v>. 

Denoting  P  ^  by  w,  and  from  the  definition  of  T(.),  this  gives 


1  T 

with  Y  ('0>  to  compute  the  action  of  the  derivative 

of  a^(u)  on  w,  we  differentiate  the  constitutive  law  (2),  first  with 
respect  to  time,  then  with  respect  to  the  velocity  u.  We  obtain 


* 

~  satisfies  the  differential  equation 

X  -  X(7u)r  -  Xr{Vu)^  +  t 


=  2pD(w)  -  X(w.V)ajj(u)  +  X^w)ajj(ti)  ♦  XOp(u)  ^7w)^ 
=  0  on  r~  (”  part  of  F  with  .  jn  <  0) . 


On  the  other  hand,  by  differentiating  the  adjoint  state  equation  (8) 
with  respect  to  time,  we  have 


(13) 


G  -  XG^u) 


+  G  ]dx 


Integrating  (13)  by  parts,  and  taking  the  incompressibility  constraint 
div  u  «  0  into  account,  (13)  yields 


n  n 


x^u);^  +  1  ]dx 


which  from  (12)  is  equivalent  to 


(14)  /  ;^.D(;^)dx 

n  ~ 


/  G.{2uD(w)  -  X(w,V)a 

n~  ~  ~  - 


Plugging  (14)  back  in  (11)  finally  gives  (10)  and  our  proof  is  complete. 


n 


With  this  theorem,  the  computation  of  J(v)  and  J'(}^  reduces  to  the 
following  sequence  of  operations. 

’ 

(i)  Solve  Pw  ■  V  in  H,  w  C  V  : 

(ii)  Compute  £p(»*ui)  by  integration  of  the  constitutive  law  (2)  ; 

(iii)  Compute  the  right-hand  side  jr  of  the  state  equation  (5); 

(iv)  Compute  the  state  vector  ^  by  solving  the  linear  problem  (5)  ; 

(V)  J(v)4  (Pj.  Px)-l<r,  2>  ; 

(vi)  Compute  the  adjoint  state  G  by  integration  of  (8) ; 

(vii)  Compute  the  right-hand  side  of  the  gradient  equation  (10)  ; 

(ix)  Compute  J'(v)  by  solving  the  linear  problem  (10). 

In  suomary,  for  the  least  squares  formulation  (5),  the  computation  of  J  and 
of  J'  requires  two  integrations  in  time  and  the  solution  of  four  linear  problems 
associated  to  the  fixed  operators  P  or  P^. 

5.  A  TYPICAL  CHOICE  FOR  THE  DESCENT  METHOD  AMD  FOR  THE  PRECONDITIONISG 
OPERATOR 


We  still  consider  the  framework  of  §3.2  and  we  now  suppose  that  V  is 
approximated  by  a  finite  dimensional  space  with  basis  (q7.)i>l,N,  to  which 
we  associate  the  matrix. 


""■^^j^-l.N.j-l.N  '  A..-/^2  u  £(jei).D(jg.)dx 


We  then  define  P  as  the  Choleski  factorization  of  A  (P^P-A,  P  lower  trian¬ 
gular),  H  as  R  and  (.,.)  as  the  canonical  scalar  product  on  . 


With  this  choice  of  P  and  H,  the  minimization  problem  (5)  can  be 
solved  by  the  standard  Polak  RibiSre  conjugate  gradient  method  given  below  : 
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given  ; 


o 

X  V 

X  jg°  •  *s  computed  in  §4  ; 


o  o 

*  «  5 


for  n>0,  with  v**  and  w**  known,  and  until  satisfied,  do  7 


X  solve  J(v  -  » 

(use  quadratic  interpolation,  for  example) 


n+‘  n  n 

XV  -  V  -  p  w 


n+ 


X  =  J'Cv***  )  as  computed  in  §4  , 

.n+*  n+1  ^  n  /  n+*  n+*  n^  ,,  n  n. 

*J1  =&  *  S.  (&  •  &  -&  >/<8  .  £  ) 


The  numerical  results  given  below  were  obtained  by  solving  (S)  with 
this  Polak-Ribiere  algorithm.  To  accelerate  convergence,  we  have  recently 
replaced  this  algorithm  by  the  one  described  in  BUCKLEY  and  LENIR  [  1983], 
which  begins  by  a  few  steps  of  a  BFGS  method  and  which  then  switches  to  a 
conjugate  gradient  method  using  the  last  BFGS  update  of  the  Hessian  as  an 
additional  preconditioning  matrin.  Again,  the  practical  implementation  of 
this  last  method  only  requires  the  knowledge  of  J  and  J’,  as  obtained  from 
§4,  and  only  involves  the  solution  of  linear  systems  which  are  associated 
to  fixed  positive  definite  symmetric  sparse  matrices  and  which  are  thus 
cheap  to  solve  even  for  N  large.  On  (S),  this  last  algorithm  gives  coBq>a'- 
rable  results  but  appears  more  robust  chan  the  original  Polak-Kibiire  me¬ 
thod. 

6.  NUMERICAL  IMPLEMENTATION 


In  the  framework  of  §3.2,  the  implementation  on  a  computer  of  the  tech¬ 
niques  of  §4  and  §5  requires  the  solution  of  two  numerical  problems  : 

(i)  what  type  of  approximation  can  be  used  for  Che  space  V  of 
velocity  fields  ? 

(ii)  for  a  given  velocity  field,  what  numerical  technique  can  be  used 
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for  the  integration  in  time  of  the  constitutive  law  and  of  the  adjoint 
state  equation? 

Those  problems  are  strongly  interconnected  since,  for  example,  the  fi¬ 
nite  element  which  is  used  deteriaines  the  aspect  of  the  computed  trajecto¬ 
ries.  D.  Malkus  [  1984]  proposes  answers  which  are  very  attractive  because 
they  respect  the  physical  structure  of  the  problem.  His  technique  decomposes 
as  follows  : 


(i)  choice  of  an  exactly  incompressible  piecewise  linear  finite  element 
(such  as  the  linear  crossed  triangle)  for  approximating  the  velocity  field; 


Hath  of  th«  doMin* 


Figure  2  -  Linear  Crossed  Triangles 
(ii)  exact  computation  of  the  trajectories  incoming  at  the  center  of 
each  finite  element  through  a  piecewise  analytical  solution  of  the  ordinary 
differential  equation 


(iii)  computation  of  the  deformation  gradient  history  by  solving  analy 
tically  the  equation 


Ft(Xt(x.T))  -  VujX^(x,T)lJ^(X^(x,T)).  F^(x,t)  -4,  ; 
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(iv)  computation  of  the  added  stresBea  by  a  Laguerre  type  numerical 

quadrature 


mj(T) 


NT 

W(T^)mj(T^)  (F^(x,t.))  . 


The  numerical  quadrature  of  (iv)  slightly  changes  the  constitutive  law 
but  respects  its  objectivity  since  the  trajectories  and  deformation  gradients 
are  exactly  computed. 


For  a  viscoelastic  fluid  with  a  more  complicated  differential  constitu¬ 
tive  law,  the  Laguerre  quadrature  is  replaced  by  a  forward  numerical  integra¬ 
tion  of  the  constitutive  law  on  the  computed  trajectories,  using  an  automatic 
time  stepping  strategy.  In  the  case  of  an  evolution  problem,  the  integration 
in  (iv)  is  only  performed  between  times  t^  and  which  is  then  much  cheaper. 

f 

In  sumnary,  if  we  use  D.  Malkus  ideas,  our  numerical  technique  for  the 
solution  of  (l)-(3)  finally  reduces  to 


1)  the  transformation  of  the  original  equations  (t)~(3)  into  an  equiva¬ 
lent  minimization  formulation  (§3). 

2)  the  solution  of  this  minimization  problem  by  a  standard  descent 
technique  (§5), 


3)  the  approximation  of  the  velocity  fields  by  linear  crossed  -  trian¬ 
gular  finite  elements  (§6), 

4)  the  construction  of  a  preconditioning  operator  by  a  Choleski  facto¬ 
rization  of  the  Stokes  operator, 

5)  the  numerical  integration  of  the  constitutive  laws  by  a  numerical 
quadrature  on  the  analytically  computed  trajectories  (§6) . 


As  described  the  method  is  expected  to  be  stable,  even  if  there  is  change 
of  type  (this  is  typical  of  a  dual  least-squares  approach),  and  accurate  (ve¬ 
ry  little  approximation  appears  in  the  integration  of  the  constitutive  laws). 
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Horeover,  although  mainly  illuatrated  in  the  upper-convected  Maxwell- 
case,  the  method  can  be  applied  to  any  constitutive  model  which  reduces 
to  a  differential  equation.  In  particular,  for  such  models,  the  gradient 
J'(.)  can  still  be  computed  as  in  §4,  simply  by  replacing  (13)  by  the 
transpose  of  the  differential  equation  defining 

Nevertheless,  the  method  appears  to  be  very  expensive,  especially 
in  the  steady  case  where  almost  all  the  computation  time  is  used  for 
the  exact  calculation  of  the  trajectories  at  a  given  velocity  field.  In 
that  respect,  the  mixed  formulation  of  §3.3  appears  more  attractive 
because  there  no  such  calculation  is  needed. 

7.  NUMERICAL  RESULTS 

We  consider  below  the  numerical  study  of  the  plane  stationary  flow 
over  a  slot  of  an  upper-convected  Maxwell  fluid.  The  domain  and  the  boun¬ 
dary  conditions  of  the  flow  are  chose  of  Fig.l,  and  its  associated  Deborah 
and  Reynolds  number  are  respectively  De>.75  and  Re>.86  1 0'** (Deborah  number  « 
product  of  the  fluid  relaxation  time  by  the  shear  rate  of  the  incoming 
fluid  at  Che  solid  wall) .  The  aspect  of  the  streamlines  inside  Che  slot 
and  of  Che  hydrostatic  pressure  profiles  are  represented  on  Fig. 3  and  4. 

The  computation  was  done  on  a  Cray  ),  with  the  finite  element  mesh 
represented  on  Fig. 2  (961  nodes),  and  using  the  Polak-Ribigre  algorithm  of 
§5.  Fifteen  iterations  were  required  for  an  execution  time  of  12  mn.  More 
than  9ST  of  this  time  was  devoted  to  the  analytical  computation  of  the  tra¬ 
jectories,  CO  be  done  twice  per  iteration.  This  indicates  that  our  stra¬ 
tegy  for  Che  computation  of  the  trajectories  should  be  revisited  in  the  sta¬ 
tionary  case,  and  that  it  may  be  better  there  to  use  Che  mixed  formulations 
of  §3.3. 


.75  ,  RO  =  .332330t-3 

a  < I  0  a  15 


Figure  3  -  Streamlines  inside  Che  slot 
(rising  flow) 


SLOT  f  r  « 1 1 1  11 

OE  =  .  75  ,  Ri  =  .  000332338 

I ( I r  « 4 1  0 1  15 

Figure  4  -  Hydrostatic  pressure  profiles 
(rising  flow) 
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